A Systematic Expansion Method in Granular Hydrodynamics 



J. WakouF] 

Instituut voor Theoretische Fysica, Universiteit Utrecht, P.O.Box 80.195, 3508 TD Utrecht, The 

Netherlands 

Abstract 

We propose a systematic expansion method which is applied to freely evolving granular 
fluids contained in sufficiently small systems. Restricting ourselves to small systems, we 
show that there exists a small parameter which characterizes a typical size of density and 
temperature inhomogeneities. A solution of the hydrodynamic equations for a fluid of 
inelastic hard spheres is expanded in the small parameter. It is shown that our method 
quantitatively describes the asymptotic state of the system, such as flow profiles, density 
and temperature inhomogeneities, and the decay law of the global temperature and the 
energy per particle. 
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1 Introduction 

Granular matter shows fluid like behavior under certain conditions |2| , such as Couette 

flow, convection flow under vibration, and flow down an inclined chute. A freely cooling 

rapid granular flow is an idealized limiting case of granular fluids in the absence of gravity 

and without any energy input. By analyzing their molecular dynamics (MD) simulations of 

a system of smooth inelastic hard discs, Goldhirsch and Zanetti |^, Q discovered that the 

freely cooling granular fluids exhibit interesting instabilities. When the system is prepared 

1 The present address is: Miyakonojo National College of Technology, 473-1 Yoshio-cho, Miyakonojo, 
Miyazaki 885-8567, Japan; E-mail: wakou@miyakonojo-nct.ac.jp. 
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in a spatially homogeneous state, it slowly develops patterns in the flow field (vortices) and 
in the density field (clusters). 

A linear stability analysis ||, |], ||, ||, |?], |8| of hydrodynamic equations for inelastic hard 
spheres (IHS) has revealed that the initial spatially homogeneous cooling state is unstable to 
the formation of vortices and clusters, and provided a good description of the initial growth 
of unstable hydrodynamic modes. At large times, however, these unstable hydrodynamic 
modes grow so large that nonlinear hydrodynamic effects, such as convection and viscous 
heating, may dominate the evolution of the system. Nonlinear effects in freely cooling 
granular fluids have been discussed by several authors ||, ||, |9], [H], (T^j . In Refs. ||, |4|, 
the importance of viscous heating effects has been pointed out and the threshold of the 
clustering instability has been given in terms of the coefficient of restitution and a typical 
length scale of a fluctuation. Nonlinear coupling between hydrodynamic modes through 
viscous heating has been studied || by means of numerical analysis of a hydrodynamic 
model and direct Monte Carlo simulation (DMCS) of the Boltzmann equation. Recently, 
strong numerical evidence has been shown [|l0| that a one-dimensional freely cooling granular 
gas asymptotically behaves as a sticky gas, and that the inviscid Burgers equation is an 
appropriate continuum theory. These three works are concerned with unstable growth 
regime in large systems, which is the physically most interesting problem of hydrodynamics 
for the freely cooling granular fluids. This regime is characterized by intrinsic patterns in 
the flow field and in the density field where the typical correlation length of these patterns 
are small compared to the system size. Unfortunately in this regime only analytical large 
time results have been obtained from a mode coupling theory (l^, 14, [D|. 

Recently some progress in the analysis of nonlinear hydrodynamic equations has been 
made for freely cooling granular fluids contained in systems that are sufficiently small, such 
that the clustering instability is suppressed [11,12]. A linear stability analysis 0, & [6|, [?|. |8| 
predicts that if the smallest wave number of a system, k$ = 2tt/L, is smaller than a critical 
value k*^, the homogeneous cooling state becomes unstable and the flow field develops a shear 
flow pattern. Several groups have observed by MD simulations ]3|, [|, ^, [L6|, 11, 17] and by 
DMCS H of IHS contained in two-dimensional systems with periodic boundaries that there 
is a regime of fco for a given coefficient of restitution where the system develops shearing 
patterns, but no strong density inhomogeneities. In this regime the flow field eventually 
reaches a stationary shearing pattern that has the following properties: (i) a stationary 
shear flow profile with a period equal to the system size in the direction perpendicular to 
the direction of the flow || f|, [f| 16, 11, |l^, [l7|], (ii) exponential decay of the energy per 
particle as a function of r, the average number of collisions suffered per particle in time 
t @, H |l|, Hi, (iii) t~ 2 -decay of the energy per particle E2l [la], (iv) inhomogeneities in 
density & ^ ||, 0, [l2|, 17] and temperature [l^] with a period that is half the period of 
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the shear flow profile; the density concentrates into regions of the peaks of the flow profile, 
and (v) saturation of the density profile ||. The properties (i)-(iii) have also been observed 
by MD simulations in five and six dimensions |l9| . In Ref. |ll| , the evolution of shear flow 
patterns in systems with k^ ~ k\_ {k$ < k^_) has been studied. In such systems, unstable 
shear modes grow so slowly that the remaining hydrodynamic modes are enslaved by the 
shear modes, and the amplitudes of the unstable shear modes remain small. On the basis 
of this consideration, they have derived amplitude equations for the unstable shear modes 
and explained the above properties in a quantitative manner. For similar small systems, in 
Ref. |12|| a Landau-Ginzburg-type equation for a nonconserved order parameter, the shear 
rate tensor, has been derived under an approximation in which inhomogeneities in density 
and temperature are neglected. It has been shown that the Landau-Ginzburg-type equation 
describes saturation of unstable shear modes due to viscous heating effects, and describes 
the properties of shear flow profile (i)-(iii) at large times in a quantitative manner in spite 
of drastic simplification achieved by the approximation. 

This paper is also concerned with the small system regime; we propose a systematic ex- 
pansion method to obtain solutions which correspond to the asymptotic state of sufficiently 
small systems. The method is applied to the stationary shear flow profile with the proper- 
ties (i)-(v), which appears in sufficiently small systems with periodic boundaries. Although 
the asymptotic state of the small systems with periodic boudaries is not a physical state, 
we concentrate in this paper only on this case, which has been studied well so far by the 
computer simulations |4|, ^, [R], [ll], 17 and by the theories [11, The method is 
based on an expansion in the Prandtl number P r of a solution of the hydrodynamic equa- 
tions for IHS. At large times in small systems, viscous heating, which is proportional to the 
kinematic viscosity v, due to growing flow field patterns becomes dominating effects that 
induce inhomogeneities in temperature || ||, ^, 11, 0. In sufficiently small systems, the 
inhomogeneities in temperature is largely reduced by thermal conduction, which is propor- 
tional to the thermal conductivity k. In these small systems, a typical size of temperature 
inhomogeneities can be determined by these two competing effects; hence it is proportional 
to the Prandtl number P r oc v j k. We will show that by expanding the solutions in the 
Prandtl number, the inhomogeneities in density and temperature field in the small system 
regime can be described; corrections to shear flow profile, and the energy and temperature 
decay laws caused by these density and temperature inhomogeneities are obtained in a sys- 
tematic manner. This method has the approximation used in Ref. |l^] that density and 
temperature inhomogeneities are neglected as its zeroth approximation. Our method can 
describe systems with finite amplitudes of unstable shear modes, to which the method in 
Ref. [ll]] can not be applied. The idea of expanding solutions of hydrodynamic equations in 
a small Prandtl number was originally developed for classical fluids by Busse (2Cj in order 
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to describe an oscillatory instability of convection rolls that appear in a fluid layer heated 
from below. 

This paper is organized as follows. In Sec. 2, we introduce the hydrodynamic equations 
for IHS and show that we can eliminate the global temperature, defined as the spatial av- 
erage of the local temperature, from the hydrodynamic equations by introducing rescaled 
hydrodynamic variables. Then, it is shown that in the small system regime some of the 
properties (i)-(v) can be associated with the existence of a stationary solution of the hydro- 
dynamic equations for the rescaled variables. In Sec. 3, the method of systematic expansion 
in the Prandtl number is presented. We first discuss the basic idea of the expansion scheme, 
then show how the expansion scheme works up to the second approximation. A quantitative 
description of the properties (i)-(v) is given in the first approximation. It presents qualita- 
tively new predictions of the behavior of d-dimensional IHS fluids. The range of validity of 
our method is discussed in the end of Sec. 3. We end with some conclusions in Sec. 4. 



2 Hydrodynamics of IHS 

2.1 Hydrodynamic equations for IHS 

We consider a system of N inelastic hard spheres, contained in a d-dimensional cubic sys- 
tem with sides L and volume V = L d . The boundary conditions are periodic boundary 
conditions. At an initial time, t = 0, the system is prepared in a homogeneous state with 
temperature To and vanishing flow field. We assume that weekly inelastic hard sphere 
systems can be described by the hydrodynamic equations supplemented by a term T that 
expresses the rate of collisional energy loss [El], [22j , 

dtn + u • Vn = — nV-u, (1) 

12 1 
ftu + u-Vu = Vp + V • On D) + V ((V ■ u) , (2) 

mn mn mn 

dtT + u-VT = _ilv-u+— V- (kVT) 
an an 

In this paper the inelasticity is always assumed to be small. The shear rate D is the sym- 
metrized dyadic, D a/3 = (V a U/3 + VpUa - \8 a pV ■ u)/2, and A : B = J2i=i E/j=i A a $Bp a . 
The transport coefficients, the shear viscosity n, the bulk viscosity £, and the heat conduc- 
tivity k, depend on the local density and temperature. 

We observe that the total energy per particle consists of a convective and internal energy 



E — E conv + Ei n t — — I ' /r 



1 2 d m 
—mnu H — nT 

2 2 



(4) 
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with rate of change, derived from Eqs. ©-(H) 



dEi n t 1 



dr 



-pV ■ u + 2r]D : D + £(V • u) 2 - |nT 



(6) 



The terms on the right hand side of the last equation represent the work done by the flow, 
the gain of Ei nt due to viscous heating and the loss of Ei nt due to collisional dissipation. 

On the basis of kinetic theory one can derive that the rate of collisional energy loss, 
r = 2jqujT, is proportional to the collision frequency uj(n,T) multiplied by the fraction 



of energy eT lost per collision [21, 22], where 70 = e/2d = (1 — a )/2d is the inelasticity 



parameter and a is the coefficient of restitution. 
2.2 Elimination of the global temperature 

In this section, we will introduce a rescaling of the hydrodynamic variables by means of the 
global granular temperature defined as 

T(t) = ^E int (t) = ly drn(r,t)T(r,t), (7) 



or equivalently the mean random (thermal) velocity vq = y2T/m. Similarly we will use 
the global density n = N/V. 

The freely cooling hard sphere fluid has no intrinsic energy scale: there is no energy 
scale in parameters which characterize interactions between particles and interactions be- 
tween a particle and a boundary. Hence, if we rescale time by t 1— > at at fixed parameters 
that characterize interactions between particles and interactions between a particle and a 
boundary, then velocities of all particles, energy, and temperature are rescaled asv^ a~ 1 v, 
E 1 — ► a~ 2 E, and T h-> a~ 2 T, respectively, but the trajectories of particles are not affected 
by this rescaling. Therefore, we can always keep one of energy scales in the system constant 
by rescaling time in an appropriate way. This fact has been utilized by Soto et al. 11] to 



enhance accuracy of their MD simulations by keeping the total energy constant. In this 
article, keeping the global temperature constant is suitable for our purpose because the col- 
lision frequency and the transport coefficients in the hydrodynamic equations for IHS are 
given as functions of temperature and density. The rescaling of the hydrodynamic variables 
given below has been used in a linear stability analysis H [5|, ||, [?], || . 

We first introduce the dimensionless time and the scaled hydrodynamic variables 

n u ~ T 

dr = udt, n=—, u = — , T = -=■, (8) 

n Vq T 

where Q = u(n, T) is the collision frequency taken at the global temperature T(t) and 
the global density n, and r is the average number of collisions, suffered by a particle in 
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time t. As a result of the rescaling, the scaled global temperature (1/iV) / drhT always 
remains constant. The collision frequency uj and the transport coefficients {r/, £, k} in the 
hydrodynamic equations will be calculated from the Enskog theory for hard spheres, in 



principle for inelastic ones 23, ^4], 25]. Moreover, fj = r/(n, T), £ = £(n, T), R = n(n,T) 
denote the corresponding Enskog transport coefficients taken at T(t) and n. Note that in 
the Enskog theory for hard spheres, Co, fj, (, R oc Vr. Then, we get a closed set of equations 
for these scaled variables, 

d T h + u • l Vn = —h i V • u, (9) 

d r ti + u-Z Vu = -^l Vp + D L ^V ■ ( v *(n)\Jf6) +— 4v f C*(«)V?V • u 

in n V J puj n \ 

- ^{d T \nf)n, (10) 

an c v n \ 

r,*(n)r ( ~~^ , 4 C CW.fcrr, ,.,2 



a n V ' a p<jj h 

- 2 7o a/ (n)f 3/2 - (d T In T) T, (11) 

where we have introduced the following notation: 

? = ^(^VT, i = C(n)\ff, "=K*(n)\/f, (12) 

77 c « 

" = oo*(n)\[f, -*L = hfp*{h)=p{h,f). (13) 
w nl 

The functions f*(n) depend on n and the coefficient of restitution a, and t/*(1) = C*(l) = 
k*(1) = u;*(l) = 1. The constants D± = vjio and -Dr = R/cppcu are the rescaled shear 
diffusivity and the heat diffusivity, respectively, and c v = d/2m and c p = c p {n) are the 
specific heat per unit mass at constant volume and pressure. The mean free path Iq is 
defined by Iq = vq/ioq. It is important to note that the constants D±, Dt, QjpCo and Iq are 
independent of T. The macroscopic equation for the global temperature T can be obtained 
from dEi n t/dt in Eq. (||) and subsequent rescaling, and yields, 

3 T lnf = i Jdr -^fa V.u+^Dj.77*(n)Vr(D:D) 

+ 3 — C (n)\/T (V • u) 2 - 2 7o u/(n)nT 3 / 2 l . (14) 
a puj 

Hence, the time derivative of In T is a functional of the scaled hydrodynamic variables n, u 
and T, and the global temperature can be completely eliminated from the equations for 
these scaled hydrodynamic variables. 
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2.3 Asymptotic state 



Equations (||)-(11) for the freely evolving IHS-fluid are supposedly describing the growing 
patterns in the density and flow fields starting from a spatially homogeneous initial state, 
and the basic objective is to describe their asymptotic evolution for small inelasticity, where 
the time scale of collisional cooling is well separated for the free time between collisions. 
Comparison of scaled field Eqs. (|9])-(|ll|) with Eq. (O) for the global temperature suggests 
that there might be two totally different scenarios leading to an asymptotic state. One 
small system scenario, where the fields equations describe a fast process with a typical 
scale t_i_ ~ L 2 /Dj_. For time r large compared to tj_ the scaled hydrodynamic variables 
reach a stationary solution and the rate of change of the global temperature, i.e. the right 
hand side of Eq. (|l4|)> approaches a nonvanishing constant. In the large system scenario, 
the field equations contain arbitrary slow hydrodynamic modes. The system, in which the 
global temperature initially decays like exp(— 2^qt) according to Haff's law p^] , gradually 
reaches an asymptotic state with inhomogeneities in the flow field, where the viscous heating 
almost entirely compensates the collisional cooling, and where the right hand side of Eq. ( |l4|) 
becomes very small. Some consequences of the second scenario have been analyzed in the 
mode coupling calculation of Ref . (l^, Q . 

Here we concentrate on a sufficiently large time scale, t > ri, in the small system 
scenario. There are some important observations that can be deduced from the assumption 
of the existence of a stable stationary solution of the scaled hydrodynamic equations. The 
results under the approximation in which inhomogeneities in density and temperature are 
neglected have already been reported in Ref. |Hq| . The assumption of the existence of a 
stable stationary solution implies that the right hand side of Eq. (|l4]) is constant, (d T lnT) = 
—27 a < 0. This leads the following important consequences: First, for large r the global 
temperature T in Eq. ([14|) decays exponentially in terms of r, i.e. 

f oc e" 27aT . (15) 



Second, the exponential decay ( |i5|) of the global temperature implies its i _2 -decay in the 
real time t. According to the definition of r given in Eq. (||), the evolution of the global 
temperature T in terms of t becomes, 

9 tT = ^rd T f = u,f (d T In T) = ^L(d T In T)T 3 / 2 , (16) 
at y ig 

where too(n, Tq) is the collision frequency at the initial time and To is the initial temperature. 
Here we have used Q oc Vt. At sufficiently large times (<9 r lnT) can be replaced by its 
stationary value — 2^f a , i.e. 

d t f = -^=2 la T*l\ (17) 
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This equation implies that T decays as t 2 at large times, which has been indeed observed 
in simulations. More precisely 

T T 



T(t) 



r 



(18) 



where t e is an integration constant. It should be noted that the prefactor of t~ 2 is a constant 
which is independent of transient behavior of the system and of the initial temperature. We 
compare this decay with the initial decay in the homogeneous cooling state (Haff's law) 

T T 



T(t)=T e 



-270T 



7o w o 



(19) 



(l+7o^) 2 

where the decay rate in r is larger because 70 > 7 a [|l^]- Third, because the collision 
frequency u is proportional to the square root of the global temperature u oc vT, the 
relation (|8]) between r and t in the asymptotic state is given by 

1 



dt 



x 



dr oc e 7aT (ir. 



(20) 



Finally, it can be shown that a stationary solution of the scaled hydrodynamic variables 
implies that the energy per particle is proportional to the global temperature: 



E 



— I dr 

N 



m 2 d ' 
— rail H — nT 
2 2 



T 



1 f , ~~2 d 

— I arnu + — 



(21) 



These results seem to be consistent with the simulation results (i)-(v) in Sec.l. Hence, 
we shall assume that for sufficiently small systems Eqs. ®-(0) have a stable stationary 
solution, and concentrate on how to obtain the stationary solution and the exponent 7 a by 
means of an systematic approximation method. 

It is important to note that the above argument is valid also for systems surrounded 
by inelastic walls as well as for systems with periodic boundaries; both systems have no 
intrinsic energy scale. 



3 Systematic expansion method 
3.1 Basic idea 

In this section we present the basic idea of a systematic expansion method to solve Eqs. 
©"(0) f° r the stationary state, i.e. with the vanishing time derivatives. In this argument 
we shall assume that viscous heating plays an essential role in inducing inhomogeneities in 
the density and temperature fields on the basis of the previous studies [||, |4], ||, 11, 12] which 



have shown importance of viscous heating effects in the formation of inhomogeneities. It 
will be shown later that the result obtained by the systematic expansion method indeed 
supports this assumption. 
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The heat locally generated by viscous heating effects is redistributed by convection and 
thermal conduction, and partially compensated by collisional dissipation. A systematic 
expansion method to be presented in this paper works when thermal conduction plays a 
dominant role in the reduction of temperature inhomogeneities through the redistribution 
and compensation of heat. The conditions necessary for this situation to be realized will be 
examined below. 

We shall first discuss a typical size of temperature and density inhomogeneities when 
effects of convection and collisional dissipation are small compared with thermal conduction 
effects, so that they may be neglected. Denoting a typical size and length scale of the scaled 
flow field are given by U and l u , respectively, the viscous heating term in Eq.(|Tl|) is given, 
in order of magnitude, by 

l D± tMy/f (d : d) ~ -D ± L- 2 U 2 (l + 0(Af,Anj) , (22) 

where AT and An are a typical size of T — 1 and a typical size of h — 1, respectively, and 
assumed to be small. We shall show later that this assumption is valid for U 2 < 1. The 
thermal conduction term in Eq.(|i~l|) is similarly given by 

£ED T Iv • («*(n)V?VT) ~ —Dt1^ 2 AT (l + 0(AT, An)) , (23) 

c v n \ / c v \ ' 

where It is a typical length scale of temperature inhomogeneities. Because the temperature 
inhomogeneities are induced by viscous heating, which is proportional to the square of Vu, 
It can be approximately related to l u by It ~ i u /2. The typical size of the inhomogeneous 
part of the scaled temperature AT is then obtained from the balance between these two 



terms (p2j) and (23), i.e. 



AT ~ ——!-lJ 2 (\ + 0(AT, An)) , (24) 
c p a V / 

where P r = D±/Dt is the Prandtl number. A typical size of the inhomogeneous part of 
the scaled density can be estimated as follows: In the asymptotic state we may assume that 
the pressure is homogeneous as a consequence of the mechanical balance. Hence the typical 
size of the inhomogeneous part of the scaled density An is related to AT by 

An = Tp B Af = p Af, (25) 

where the constant (3b is the bulk expansion coefficient; its definition is given in Appendix 
A. It is shown in Appendix A that for inelastic hard sphere fluids (3q = T(3b is a constant 
independent of T, and in two and three dimensions (3q is less than 1, i.e., Ah < AT. 

Values of (c v /c p )(P r /d) for elastic hard spheres of various densities are listed in Table 1. 
As shown in Table 1, the factor {c v /c p ){P r /d) is a small quantity and hence AT <C 1 
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Pr 


Cp/ c v 


(c v /c p )(P r /d) 


dilute limit (2D) 


0.5 


2 


0.12 


dilute limit (3D) 


0.67 


5/3 


0.13 


Enskog formula (2D, = 0.1) 


0.46 


2.0 


0.11 


Enskog formula (3D, <p = 0.1) 


0.60 


1.7 


0.12 


Enskog formula (2D, <fi = 0.4) 


0.44 


2.2 


0.11 


Enskog formula (3D, = 0.4) 


0.89 


2.4 


0.12 



Table 1: Values of P r , c p /c v and (c v /c p )(P r /d) given by the Boltzmann theory (dilute limit) 
and by the Enskog theory (Enskog formula) for elastic hard spheres in two dimensions (2D) 
and in three dimensions (3D). Two values of the packing fraction, <f> = 0.1 and 0.4, are used 
for the Enskog formula. 



if U 2 ~ 1. This result justifies the approximation used in Ref 12], in which density and 
temperature inhomogeneities are neglected, for a finite size of shear flow profile with U 2 < 1 
in the asymptotic state. 

The fact that the inhomogeneous part of the scaled temperature and density is small, 
and is proportional to the Prandtl number motivates to expand the stationary solutions 
in the Prandtl number, expecting that it will eventually lead an expansion in the small 
parameter (c v /c p )(P r /d), which characterizes the small magnitude of density and temper- 
ature inhomogeneities when U 2 < 1. It will be shown that up to a part of the second 
approximation studied in this paper, this is indeed the case. 

In Sec. 4, we will discuss the range of validity of this method and give a system size 
above which our analysis breaks down. It is noteworthy that temperature profile produced 
by the combination of viscous hating effects and thermal conduction effects appears also in 
flows of classical fluids such as the compressible Couette flow for a fluid confined between 
isothermal or adiabatic walls [26]. 

Finally, let us consider the conditions under which thermal conduction may be regarded 
as dominant in comparison with convection and collisional dissipation. The condition that 
the convection term in Eq. (pT|), Iqu ■ VT ~ IqUI^AT, is negligible compared with the 
thermal conduction term yields, 

_1 Af / c„ P r \ d 



l()Ulrp 



^D T l^ 2 AT 



c p d 



Re <C 1. 



(26) 



Here R e is the Reynolds number defined by R e = Ul u /u, where U = vqU is a typical size 
of the flow field (not scaled). Smallness of the parameter (c v /c p )(P r /d) suggests that this 
condition is fulfilled if 

d 

— j 

2 



^Re £ 1. 



(27) 
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The inhomogeneous part of the collisional dissipation term, which contributes to reduction 
of inhomogeneities in the temperature field, is in order of magnitude given by 70AT. The 
condition that this term is small compared with the thermal conduction term yields, 

7Q A ^ _ [ 2l^A d ^° 1 (oa\ 

Since (c v /c p )(P r /d) <C 1, this condition is fulfilled if 

< 1. (29) 

For the stationary solutions with the properties (i)-(v), which are to be studied in this 
paper, however, the property (iv) and the assumption of homogeneous pressure suggest that 
temperature are inhomogeneous only in directions perpendicular to the direction of flow; 
hence the convection term identically vanishes. Therefore the condition (p7|) is not required 
for this type of solutions. Besides, it will be shown in Sec. 3. 4 that for solutions with the 
properties (i)-(v), the condition ( |2^ ) is fulfilled if U 2 ~ 1. 

On the basis of the consideration denoted above, we expand the stationary solution in 
the Prandtl number as follows: 

n = 1 + en (1) + eV 2) H , 

u = u<°> + eu« + e 2 u^ + . . . , 

f = 1 + eT (1) + e 2 T (2) + • • • . (30) 

Here we have introduced a parameter e that is of 0(P r ). 

We require that the conservation of the total number and the total momentum, N = 
J dm and = / drnu, is satisfied in each order in e. This gives rise to the following 
conditions: 



= / drnW for I = 1, 2, - - (31) 
= [dru®, 0= /drffiW+nWuW), (32) 

In addition, rewriting the definition of the global temperature as 1 = (1/V) / drhT and 
requiring that this relation is satisfied in each order in e, we get the following conditions for 
the scaled temperature: 

= [drfW, = Jdr(^fV+f^), .... (33) 



We substitute ( |30| ) into the hydrodynamic equations (p|)-(pd|) with the vanishing time 
derivatives 

u-Vn = -nV-u, (34) 
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nu • IgVu 



nu • l VT 



Vp+2D±V ■ n*(n 



TD) + -=V [C(n)yTV u 



- (<9 T lnT) nu, 

a c„ 



(35) 



K*(n)\JTVT 



+ 



d 



D±r]*(h] 



/' f D : D ) + - — C{n)\ff(V ■ u) 2 
a puo 



2-f uj*(fi)fif 3/2 - (d T InT) nf. 



(36) 



Here both sides of Eqs. (10) and (|TT| ) are multiplied by n for later convenience. We put e _1 
in front of the thermal conduction term, because P~ l appears there by rescaling collision 
time r by L 2 /Dj_ and by rescaling space r by L. We omit this procedure since this rescaling 
of variables by constants does not change the final result. 

Collecting terms of the same order, we obtain a set of equations that determine a 
stationary solution of the scaled hydrodynamic variables in each order in e. The smallness 
parameter e is set equal to 1 in the end of calculation. 



3.2 Zeroth approximation 

The terms of 0(e _1 ) gives an identity = 0. The zeroth approximation, determined 
by the terms of O(e ) 



i.e. 







u(°) • l Vu(°) 



D ± V 2 u(°)-i(5 T lnr) (0) uW, 



where 



(drlnT 



(0) 



1 

V 



— \ dr 



- d D L (Vu(°))t:Vu(°)-27o 



(37) 
(38) 

(39) 



(|38| ) has only the trivial solution u(°) 
of the system and k*. 



where (Vu^)^ = dpu^ ■ We show in Appendix B that if k$ > kj_, the nonlinear equation 

0. Here ko = 2tt/L is the smallest wave number 
\rjo/D± is the critical wave number of the instability in shear 
modes, predicted in a linear stability analysis (3|, |], ||, ^, |7], B. In this case all higher order 
corrections vanish, suggesting that the system is at large times in a homogeneous cooling 
state that follows Haff's law |J7|, (<9 r lnT) = (<9 T In T) ^ = —2jq, independently of the 
initial state of the system. 

Hereafter, we shall consider only the non-trivial case ko < k*^. If ko < k*^_, there is a 
large set of solutions |12j that satisfy Eqs. (p?|)-(|39|). In this paper we do not pursue to find 
all possible solutions that satisfy Eqs. (|37|)-(|39|). We show that there is indeed a solution 
that corresponds to the flow profile observed in simulations [||, ||, [6|, [16], 11, |l^, [HJ with a 
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period of the system size L. An intuitive explanation of the flow profile with the period L 
is that this is the mode which shows the slowest hydrodynamic decay 0, |4| . Assuming the 
flow profile on the basis of the observations, we concentrate on a quantitative description 
of the amplitude of the flow profile, which eventually leads us to a theoretical prediction of 
the amplitude of inhomogeneous temperature and density profile. 

One can show that the solutions of Eqs. (|37|)-(39) with a period of the system size 



L have to be a flow that is inhomogeneous in directions which are perpendicular to the 
direction of the flow [see Appendix C], i.e. 

d! d 

u {0 \r) = J2 A a pe/3cos(kor a + 6 al3 ), (40) 

a=l p=d'+l 

where we have chosen a = 1, 2, • • • , dl (dl > 1) as to be directions in which u(°) is inhomoge- 
neous, and /3 = d! + 1, • • • , d as the remaining directions, is a unit vector in the direction 
(3, and a p is an arbitrary phase factor. The amplitudes A a p are real numbers that satisfy 
the following conditions. First, from Eq. (j37|), we get A aa = 0. Second, if ko < k]_, Eq. 
give rise to a relation |l2| : 



4 = E t ^ = db °- D / l) - (41) 

The proof of the results (fl(]) and (fll|) are given in Appendix C. From Eqs. ( f40|) and ( fill) 
we obtain a typical size of the scaled flow field U = [(1/V) J dr\u^ (r)] 2 ] 1 / 2 = A /y/2. 

The achievement of one of the stationary solutions can be viewed as a process of spon- 
taneous symmetry breaking fli~2| . Good quantitative agreement about the amplitude Aq 
between the theoretical prediction ( f4l"D and the result of MD simulation has been shown in 



Ref. p|. 



Substituting the solution (||) into Eq. (||), we get 

(d r lnf) (0) = -2D ± kl (42) 

3.3 First approximation 

First order correction can be obtained from the terms of 0(e°), 

v2 ~ (1) = _ Cjl a r (0) )t . v -(o)_ i f d r(va( o) )t : va( )l, (43) 

c p d l V J J 

where we have used a relation Vu' ' : Vu' ' = 0, which is obtained from Eqs. ( |37|) and 
The relation (E^) manifests that the temperature inhomogeneities are excited by 
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the inhomogeneous part of the viscous heating term. Substituting the expression © and 
solving the Laplace equation under the condition (33), we get 

T (1) (r) = -f^lX 2 cos [2 (kor a + 9' a )], (44) 

P a=l 



where A' a and 6' a are defined by A'£ cos [2 (k r a + #„)] = J2$=d'+i ^a(3 cos [ 2 (^o r a + #a/3)]- 
It can be shown that if d > 3, there exist solutions with the constants ^4 a ^ and 9 a p with 
vanishing However, in the following analysis we will disregard these solutions because, 
to the best of our knowledge, the asymptotic states corresponding to these solutions have 
not been reported. Then, we will chose B' a so that A'£ = J2p=d'+i cos P ~ @' a )] 
becomes positive. 

The first approximation of the scaled density and flow field can be obtained from the 
following relation given by the terms of 0(e 1 ): 



(45) 



u(°)-ZoVuW+uW-ZoVuW 



lo 
2 



|Q n(D + if (Vu(°) + (Vu^)t 



+ D ± V 

- i(d r mf) (0) (n«u(°)+uM) - i («9 T lnT) (1) a 



va( 1 ) + (Vu«)t-±[/(v-iiW; 



(0) 



(46) 



2 v 7 V 7 2 

where (• • -)o represents that (• • •) is evaluated in the homogeneous state, i.e. n = n and 
T = f ; (d T lnf) (0) is given by Eq. @ and 



(5 T lnf) 



(i) 



V 



r/r 



<9n /„ 



+ -Di(Vu(°))t:V# 
d 



(47) 



These equations can be solved analytically for ra'W and u'W that are spatially averaged 
over the directions j3 (/3 = d' + 1, • ■ ■ , d), in which u(°) is homogeneous: 

- {1) = jh*l A ^ = j^j A (48) 

J /3=d'+l /3=d'+l 

By taking the spatial average over the directions /3 (/? = d' + 1, • • • , d), of both sides of Eqs. 
and (|4"6|), we get a closed set of equations for n W and u W: 



(49) 



«=i 
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lo 
2 



dnj 



+ D ± 



/ dr]* 
\ dn 



\dfJo 

V 'n'W + Vf( 1 ))-V'« ( i )) + 



(50) 



^) o n (1) + ^(^V' 2 ^ 0) 



+ D ± V 2 u'^ + D ± k 2 u'W + D ± ifegn'«fiJ ,) - ~ (d T lnf uf, (51) 
where a = I, ■ ■ ■ , d' and /3 = d' + 1, • • • , d, and V' has been defined by V' = (di, d%, ■ ■ ■ , 8^)- 



Here Eq. (42) has been used and 

d' 

(d T \nff 



^/^*4 D 4((f) " (1,+ H (v ' i<0,,,:V " < ° , 



+ 2(V'ii( ))t:V'u« 



The equation (|50|) can be solved without any difficulty. Combination of Eqs. (|49|) and (pQ 
yields 



(52) 



<9n/ 



(53) 



Because of the periodic boundary conditions and the conditions ( |3l| ) and (p3|), we obtain 

(§) 



n'Wfr) 



(H), 



^Vor( 1 )(r) = -^ T«(r). 



I (§E 

T \dn 



(54) 



This relation is consistent with the pressure balance achieved in the asymptotic state. Sub- 
stituting Eqs. fl44j) and fl54| ) into Eq. (|50|), we get 

= D ± V' 2 ul (1) + D ± fcgn^ 1 ) . (55) 

This relation suggests in terms of the Fourier modes tt^ k = if |k| ^ ko; hence u a 
(a = 1, ■ ■ ■ , d') is given by 



*4 (1) ( r ) = H Cq'q COs(/c ^a' + C*'a), 



(56) 



a'=l 



where the amplitudes C a > a are real numbers that satisfy C aa = because of Eq. (fH), and 
the phase factors ( a > a are also real. 

Since T^-\ and uj 1 ' (a = 1, • • • , d') are now known, Eq. (|5l|) is an inhomogeneous 



linear equation for Un (/3 = d' + 1, • ■ • , d). We will show in Appendix D that the solubility 



(57) 



conditions of Eq. ( |5l|) is satisfied if and only if 



t or 0' a + ir, £ Al p = X 2 = ^Al 



P=d'+1 
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for a = 1, ■ • • , d! and (3 = df + 1, • • • , d. The freedom of the choice 9 a p = 9' a or 6' a + tt can 
be absorbed into the definition of A a p\ hence Eqs. ( f40|) and (44) can be rewritten as 



u(°)(r) = Yi, Yl A a pe cos(k Q r a + 9' a ), 

a=l (3=d'+l 

rfi d' 



- J2 cos [ 2 ( + O'a 



d 



(58) 
(59) 



a=l 



where T = — {c v /c v ){P r /2d)A 2 ] . According to the relation (54), n W i s given by 

* d' 



(60) 



where n = — j3$T. Once the solubility conditions are satisfied, calculation of the solution of 
Eq. (|5l"|) which satisfies the condition (32) is straightforward. The final result is presented 
in Appendix E. 

Finally, substituting Eqs. (f|), flU), © and flA~35l ) into (||), we obtain the first order 
correction to (<9 T lnT), 



(d T In T") 



(i) 



^4 



1 + 



9r/* 
<9n 



h+ -T 
2 



(61) 



Hence, from Eqs. ( f42|) and (|6lD, the exponent 27 a in Eq. ( |i~5| ) is given by 



2 7 a 



(<9 r lnT) (0) + (<9 T lnT) 

V 



(1) 



(62) 



Notice that the first order correction to 2^ a depends explicitly on the coefficient of restitution 
a through 70. The ratio of the energy per particle to the global temperature is independent 
of time and can be calculated from Eq. ( ^T|) to order e: 



E 

drp 

2 1 



2 1 

dV 

7o 
D ± k 2 



J dr (xi^ 2 + 2u(°) • u« + ft^u^ 2 ) + 1 



1 + ^ P r (70 ~ D ± k 



D±k { 



(63) 



where Eqs. @), © and (|A3^ ) have b een used. 

Calculation in the higher order approximation is rather involved. In Appendix F we 
discuss the second order correction to the scaled temperature, T^ 2 \ in order to illustrate 
how excitation of modes with larger wave numbers can be described in the higher order 
approximations . 
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3.4 The range of validity of the expansion method 

A restriction to the range of validity of the expansion method is given by the condi- 
tion that the amplitude of the temperature inhomogeneities \T\ = {c v /c p ){P r /2d)A^ = 
(c v /c p )(P r /d)U 2 , which is greater than the amplitude of the density inhomogeneities |n| in 
two and three dimensions, has to be much smaller than 1. In Sec. 3.1 we have seen that 
(c v /c p )(P r /d) is a small parameter; therefore the above condition is fulfilled if 

~ 2= Aj = d l0 -D ± k 2 

2 2 D ± k 2 V ; 

This condition can be rewritten as 

70 < (l + ^D_Lfcg, (65) 

and hence tells how small the system has to be in order that our systematic expansion 
method describes well the asymptotic state. 

We have also seen in Sec. 3.1 that the inhomogeneous part of the collisional dissipation 
term, which contributes to reduction of inhomogeneities in the temperature field, is neg- 
ligible if the condition (^) is fulfilled. From the condition (|65|), however, we see that if 
U 2 < 1 then the condition ( p9| ) is fulfilled. This result is consistent with the fact that the 
inhomogeneous part of the collisional dissipation term appears in the second order correc- 
tion to the scaled temperature, as shown in Appendix F. Therefore the range of validity of 
our systematic expansion method applied to a solution with the properties (i)-(v) is given 
only by the condition U 2 < 1. 



4 Conclusion 

In this paper, a systematic expansion method in the Prandtl number has been developed 
in order to describe the asymptotic state of freely cooling granular fluids contained in suffi- 
ciently small systems. Although the Prandtl number of IHS fluids itself is not a very small 
parameter, it eventually leads an expansion in the small parameter (c v /c p )(P r /d), which 
characterizes the size of temperature and density inhomogeneities if a typical size of the 
scaled flow field U = A /V2 satisfies U 2 < 1. The expansion in the Prandtl number yields 
a set of inhomogeneous linear equations. The zeroth approximation is a homogeneous state 
where only the scaled flow field is allowed to have finite amplitude. The resulting equa- 
tions in the zeroth approximation are in agreement with those equations for the stationary 



state which are given in Ref. [12| with an approximation in which density and tempera- 
ture inhomogeneities are neglected. It has been shown that inhomogeneities in density and 
temperature are systematically described in the higher order approximation; corrections to 
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macroscopic quantities such as the energy per particle and the global temperature caused 
by these inhomogeneities are estimated up to the first approximation. 

Our theory describes sufficiently small systems where thermal conduction plays a more 
important role than convection and collisional dissipation. As a consequence, the temper- 
ature inhomogeneities is determined mainly by viscous heating induced by the flow profile, 
and thermal conduction. Density profile is then determined in such a way that pressure av- 
eraged over the directions of the flow in the zeroth approximation remains homogeneous. If 
the size of system is so large that thermal conduction as a diffusion process is less important 
than collisional dissipation for large scale fluctuations of temperature, a nonlinear theory 
by Goldhirsch et al. [||, Q is more appropriate to describe temperature profile determined 
mainly by viscous heating and collisional dissipation; their theory describes also pressure 
imbalance that eventually leads a clustering process. Apparently the clustering process can 
not be described by our method because the condition of small inhomogeneities in density 
and temperature is necessary for the systematic expansion method to work well. 

The condition that the amplitude of temperature inhomogeneities \T\ = (c v /c p )(P r /2d)A 2 l 
(c v /c p )(P r /d)U 2 has to be much smaller than 1 yields a restriction on the system size L for a 
given packing fraction and the coefficient of restitution a of IHS. Because (c v /c p )(P r /d) <C 1 
in two and three dimensions [see Table 1], our theory is expected to work well for U 2 < 1; 
this is the case that can not be described by those amplitude equations in Ref. Jll[] which are 
derived assuming £/ 2 <C 1. One can show that our results expanded in the small parameter 
in Ref. |11], U 2 = Aq/2, coincide to order A\P 2 with the results in Ref. expanded in 
the Prandtl number P r , except for the indefinite constants in our theory, which can not be 
fixed in the first approximation. 

In this paper, stability of the stationary solutions has not been discussed. We have 
assumed the flow profiles that are observed in computer simulations (the property (i) in 
Sec.l) and concentrated on the results that can be deduced from this assumption, such 
as the amplitude of the scaled flow field, the amplitude of the density and temperature 
inhomogeneities, and the decay exponent of the global temperature and the energy per 
particle. The stability analysis may eliminate solutions that are not observed in computer 
simulations and therefore not discussed in this paper. 

In this paper the systematic expansion method has been applied to freely evolving granu- 
lar fluids contained in small systems with periodic boundary conditions; these systems have 
been studied in the earlier works by computer simulations and by theories. The asymptotic 
state of these systems is unphysical because it is dominated by the finite size effect due to 
periodic boundary conditions. Our method of systematic expansion, however, is expected 
to be useful to study other practical problems which occur in finite size systems with phys- 
ical boundary conditions, such as Couette flow and flow under vibration, if the dimensions 
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of the system are sufficiently small such that viscous heating and thermal conduction play 
dominant roles in formation of density and temperature inhomogeneities. 
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A Bulk expansion coefficient for inelastic hard spheres 

In this appendix, we present the expression for the bulk expansion coefficient (5b for d- 
dimensional IHS fluids and show that (3q = T(3b < 1 in two and three dimensions. The 
bulk expansion coefficient (3b is defined by 



1 dn 

" B = 

n 61 



1 


dp 




dT 


n 


n 


dp 






dn 


T 



(A.l) 



The expression for the local pressure in Enskog's kinetic theory has been given in Refs. pl| , 
22] for three-dimensional IHS fluids. A generalization to d dimensions has been given in 



Ref. @: 



V 



nT 



i + a d n d 

1 H rna — - 

2 A 2d 



(A.2) 



where Q d = 2vr d / 2 / r (^/ 2 ) is the surface of a ti-dimensional unit sphere, a is the diameter 
of the hard sphere, and x is the pair correlation function for two particles in contact. In 
the Enskog approximation x has the same dependence on the local density at the point 
of contact as the equilibrium pair correlation function at contact has on the equilibrium 



density. Substituting the expression ( A.2 ) into Eq. (AT), we get 

1 



Pi 



("^ x dn 



2d 



From this expression, we see that (3q = T(3b is independent of T, and if 



^(Xn)>0, 



(A.3) 



(A.4) 



then (3q < 1. The condition (| 
Carnahan-Starling approximation 



X given by the Verlet-Levesque approximation |p0 |, x 
nft d (a/2) d jd is the packing fraction in d dimensions. 



is satisfied in three dimensions for x given by the 
| , X = (2 — 0)/2(l — 0) 3 , and in two dimensions for 

) 2 , where <h = 



(1 - 70/16)/(l 
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B Homogeneous cooling state 



In this appendix, we show that if k > k^, Eq. ( |38| ) has only the trivial solution u^ ) = 0. 
We first take the inner product of u^ ) with Eq. (|3^) and integrate both sides over the 
whole space: 

J dru® • (u(°) • Z Vu(°)) = D ± J dru^ ■ V 2 u^ - ± (d T lnT) (0) J dr(u^) 2 (A.5) 

Using the relation u(°) • (u(°) • Vu' ') = u(°) • V(u(°)) 2 /2 and performing an integration by 
parts, it can be shown that the left hand side of Eq. ( A.5|) vanishes as a consequence of 



Eq. ( p7\) and the periodic boundary conditions. Rewriting the right hand side in terms of 
Fourier modes defined as = / dre~* kr a(r), we get 

= E (70 - - 2 -D ± ip?|u£>| 2 j (A.6) 

We notice that if ko > k*±_ = \J 70 /D± , the inside of the parenthesis on the right hand side 
of Eq. ( |A.6[ ) is negative for any k. Therefore, the equality is satisfied if and only if u*- -* = 0. 

C Flow profile given in (|30| ) 

The flow field that consists only of the Fourier modes with the smallest wave number 



ko is written in general as \12\ 



i (0) (r) = E E A «^P cos ( + Bap), (A- 7) 

a=l/3=l 



where the coefficients A a p and a p are arbitrary real constants. From Eq. (|37|), we get 



For the flow field given by Eq. ( |A.7| ), both sides of Eq. ( |38[ ) have to vanish independently 
of each other, because the right hand side of Eq. (38) consists only of the Fourier modes 



with the smallest wave number ko, while the left hand side of Eq. (38) does not contain the 
smallest wave number modes. We first show that the condition of the vanishing left hand 
side of Eq. (38) leads to the flow field (f40|). The condition of the vanishing left hand side 



can be written in terms of the Fourier modes, 

= E (A.8) 

ki 

for any k. Because = if |k| ^ ko, the non-trivial conditions are given only for the 
wave numbers k that satisfy |k=pko Q | = ko (a = 1,2, - ■ ■ ,d), or, equivalently, only for 
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k = ko a ± ko/j (a = 1, 2, • • • , d, j3 = 1, 2, • • • , d). Here ko Q = koe a and e a is a unit vector in 
the direction a. For a given k = ko a + ko/j, the condition ( |A.8| ) is written as 








u 



(0) 



,(o) 



0/3 u k„ a + u k„« " k 0a u 



r,(0) 



u 



<();3 

i(0)* 



(0) 

ko a 



i(0> 



■ k 0/3 U knfl + U k nfl • k 0« u k ( , 



(0) 



(A.9) 
(A.10) 



*0/3 ' K 0,3 

Substituting the Fourier transform of Eq. ( |A.7| ), Uk. 0a = (V/2) Sa=i &pB<xP> where i? Q( g 
j4 0( g exp(i6' cej a), into the above condition ( [A.IOD , we obtain 



— B a pB^ + Bp a B ai , 

= BapBp^ 



Bj3 a B a ^y, 



(A.11) 
(A.12) 



for any set of (a, (3, 7). 

This conditions imply the following properties of the complex matrix B a p. First, con- 
sidering the case when 7 = a in Eq. ( A.ll| ), we get 



= BafiBfia, 



(A.13) 



for any a and (3. Here we have used property that B aa = 0. Second, suppose B a p / 0. 
Then, Bp a = from Eq. ( A.13j ) . Then, the condition (A. 11) is reduced to = B a pBp^. 
Because B a p ^ 0, we get B^ = for any 7. Then suppose A a p ^ for a set of (a, (3) = 

mi^m) 8-nd A a p — otherwise. Let us label the directions in the 
subspace spanned by fj,%, • • • , fj, m ) as (1, 2, • • • , d'), where d' is the number of dimensions 
of the subspace. Note that d' < m because some of [n may be the same direction. Then 
one finds that any can be equal to none of (1, 2, • • • , d'), because B, 



Vi7 



A 



for 



any 7. The condition ( |A.12|) gives no additional conditions. Hence, we may conclude 
that the direction of the flow, which is in the subspace spanned by {u\, 1/2, • • • , u m ) has to 
be perpendicular to the directions (1, 2, ■•■ , d') in which the flow shows inhomogeneities. 
Therefore the flow can be in general written as Eq. (40). 



Then it is straightforward to get the condition (41) by substituting the expression (fjQ| 
into the condition of the vanishing right hand side of Eq. (j38[) . 



D Solubility conditions of Eq. ( |5T1 ) 



'(11 

Equation (51) is the inhomogeneous linear equation with respect to Uo and rewritten as 

Cuf ] =fp, (A.14) 



where 



£ u 



'(i) 



D ± v' 2 uT ] + D ± ktu 



l >3 



^/n^- (o) - ,(i) 

a=l 



~(0) 

u 



(A.15) 
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and 



fp = u'W.l Vnf-D ± kin^ 



)~(0) 

>u- 



D± 
2 



dn 



d 



1 /■ d ' 



3r/* 



1^ 



™ / o 



'(i) + ^T( 1 ))(V / a(°))t:V / l i( ) 



r,(0) 



(A.16) 



where u(°), f (1 \ nW, and fiW are given by Eqs. @,©,(§3),(|05|) and Let us de- 
fine a Hilbert space |a) to represent any function ap(r) ((3 = d'+l, ■ ■ ■ , d, r = (n, V2, ■ ■ ■ , r^')) 
by 



a /3( r ) = (r,/3|a) • 
We define the scalar product of two vectors | a) and | b) by 



^ „ d! d 
(a\b) = jj /J! dr a ^ a^r^r 



(A.17) 



(A.18) 



With this definition, the inhomogeneous linear equation ( A, 14 ) can be written as 

£\u'W) = \f), (A.19) 

and the linear operator C is Hermitian: 

(a\Cb) = (b\Ca) . (A.20) 

The inhomogeneous linear equation ( |A.19[ ) has a nontrivial solution if and only if |/) is 
orthogonal to the solutions \4> l ) of the homogeneous problem C \ <j> 1 ) = |3l|], i.e. 



V) = o. 



(A.21) 



One can show that the linearly independent solutions \4> l ) of the homogeneous problem 
are given by 

d> 

<^(r) = (r,(3\cp l c ) = Y,c l ap co S (k r a + e af3 ), (A.22) 



a=l 
d! 



> s p 



a/3), 



(A.23) 



a=l 



where the coefficients c a/3 (I = 1, 2, • • • , d'(d — d') — 1) and s™p (in = 1, 2, • • • , d'(d — d')) are 

-aP 

d' d 



real numbers. The coefficients c l 3 satisfy the condition 



E E * 

a=l/3=d'+l 



(A.24) 
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We also require the orthogonality property 

d' d 

2 
1 



i\€, 



^E E C af3 C U = 5 «'> 
a=l/3=d'+l 

d' d 



I I Jjn\ 



E E 



m m 



a=l 



0. 



The solubility conditions (|A.21|) , (<#.|/) = and = 0, yield 

d' d 

E A^A a/3 c l at3 cos[2(9' Q -e a/3 )] = 0, 

a=l/3=d'+l 



(A.25) 

(A.26) 
(A.27) 

(A.28) 



and 



£ £ 4, 2 A^sin[2(^-M] = 0, 
a=l /3=d'+l 



(A.29) 



respectively. Here we have assumed [(1 + {drf /dh)o)Po — 1/2] ^ 0, which is satisfied in 
general. In deriving Eq. ( A.29Q we have used the fact that = (a = 1, • • • d!) if |k| ^ ko 
as shown in Eq. (|56|), 

Because the coefficient sig is arbitrary, let us choose it as 



4/3 oc A a psm[2(6' -6 a p)]. 
Substituting this expression intoEq. ( |Q§D , we get 

2 £ < 2 ^sin 2 [2(^-M] = 0. 
o=l/3=d'+l 



(A.30) 



(A.31) 



Because > by definition [see note below Eq. (|44|)1, the condition ( [A. 31 ) suggests that 
if A a p ^ 0, then (6 a p — 0' a ) is either 0, 7r/2,7r, or 37r/2. We note that if A a p = then the 
phase factor 6 a p does not appear in our theory [see, Eq. (09)]. If the condition ( A.31| ) is 
satisfied, the conditions ( A.2S| ) for s™p (m = 2, 3, • • • , d'(d — d')) are automatically satisfied. 

The condition ( A.28 ) has to be satisfied for a set of constants cLg, which satisfies Eqs. 
( |A.24| ) and ( |A.25[ ). This implies that A'£ A a/3 cos[2(6' a - 9 a p)\ = const, x A a p, i.e. 

A'a cos[2(9' a - Oap)] = const. = Z, (A. 32) 

where the constant Z is independent of a and (3. The condition that the left hand side of 
Eq. QA.32 ) is independent of j3 suggests that (0 a p — 0' a ) is either or tt, or (6 a p — 9' a ) is 
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either tt/2 or 3n/2. Because of the definition of 9' a [see note below Eq. flUf)], however, only 
the former choice 



~>a.p 



9' = or vr 



(A.33) 



is allowed. Then, from Eqs. ([A321) and flA~33p we get A'£ = J2p=d'+l A lp = z - Finally, 
the constant Z is determined by the condition (|4l|), 



(A.34) 



E First order correction u 



'(i) 

/8 



First order correction Ug can be obtained by solving Eq. (|5ll) under the solubility condi- 



tions (57) and the condition (32): 
d' 



ig^(r) = Z^l ^ A a/3 COs[3(/c ? , a + O] 
a=l 
61 d< 

+ U *Y,Y, A «'P { cos [ fc o(2^a - r a >) + 2^ - + cos [k (2r a + r a ,) + 26' a + 6' Q ,] } 



a=l ot'=l 

I d ' 

+ ttV E A c*pUa X) (r) sin(/c r a + 0' a ) 
D±k a=1 

d' d' 
where the coefficients E a p are real numbers that satisfy the condition 



(A.35) 



»=l/?=d'+l 



(A.36) 



and the coefficients F a p are arbitrary real numbers. Here the constants Ws are defined by 



u 2 

U 3 



1 



lQd' 
1 

_d_ 



\oa)J 2 
«V 2 



(A.37) 
(A.38) 
(A.39) 



It should be noted that all of the W's are proportional to the small parameter (c v /c p )(P r /d). 
The results ( |56[ ) and ( |A.35 ) show that u'^ 1 ) has indefinite constants that can not be fixed 
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in the first approximation. We expect that they are fixed by the solubility conditions of the 
inhomogeneous linear equations in the higher order approximations. As shown in Eqs. ( |62| ) 
and (|63|), however, the expressions for macroscopic quantities such as the global temperature 
and the energy per particle up to the first approximation do not contain these indefinite 
constants. 



F Second order correction 



The purpose of this appendix is to illustrate how excitation of modes with larger wave 
numbers can be described in the higher order approximations. 

The spatial average over the directions (3 = df + 1, • • • ,d of both sides of the 0(e 1 ) 



equation obtained from Eq. (|36|) yields 



PD T V 2 f'W = -^D T V 



d 
270 



dn* \ 
dh J 



n'O) + If W ) V'f W 
dhj 2 1 



h'U + If w 



(V'u(°))t : V'uW - ^(V'u^jt : V'u'« 



(^gi^^W + 3f(l)j + {drln ff) (n'd) +f(D) + ( 5r lnf) (1) , (A.40) 



where T is that second order correction to the scaled temperature which is spatially 
averaged over the directions (3 = d' + 1, • • • , d. This is again the inhomogeneous linear 
equation (the Laplace equation) for f '( 2 ) and it can be shown by substituting the result of 
the first approximation that the solubility conditions of this equation are satisfied. Because 
the general expression for f '( 2 ) is rather lengthy and complicated, we concentrate here only 
on the case that the zeroth approximation u(°) is inhomogeneous only in one direction, say 
a = 1, i.e. d = 1; this is always true in two dimensions if k$ < k*^_. Then u 1 = and for 
P = 2,3,---,d, 



^ (n) = UtAtp cos[3(A;or 1 + 0[)] + E lf3 cos( Vi + 0[) + F w sm(k n + 0[), (A.41) 



where the coefficients E\p satisfy J2p=2 A.\$E\p = U3. Substituting the expressions ( p8| 
(|59|), ( |60| ) and ( A.41| ) into Eq. ( A.40| ), and using the condition (^), we eventually get 



f' (2) (r!) = Ticos[4(A;ori + (9i)] + T 2 cos[2{k ori + 6[)} 
+ T 3 sin[2(fc r 1 +^)]+r 4 , 



(A.42) 



where 



T 
32 



dh 



on ) y 



13, 



(A.43) 
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* - »'+'(£)„«(£)>*f* 
+ ^[H£)„ +s (£)> + §*]}- ^ 

d 

T 3 = 2TAfY, A W F ^ ( A - 45 ) 

0=2 

T 4 = -™ (A.46) 

Here has been obtained from Eq. (|33|). It should be noted that 71, T2 and 74 are 
proportional to the square of the small parameter (c v /c p )(P r /d); T3 contains the indefinite 
constants Fip, which are expected to be fixed in the higher order approximations. Gener- 
alization of this result to the case of an arbitrary given by Eq. ( ]58| ) is straightforward. 

In order to make analytical calculation executable in the even higher order approxima- 
tions, it would be necessary to introduce a physically sensible assumption that the pressure 
is homogeneous in the asymptotic state, which is in accordance with the result (|54|). Con- 
sequences of this assumption, however, shall not be pursued in this paper. 
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